%%%%%%%%%%%%%%%
% Housekeeping
%%%%%%%%%%%%%%%

clc;
clear                                   all;
close                                   all;
beep                                    off;
rng.Seed                                = 110887;
set(groot,'defaultAxesFontName', 'Garamond')
set(groot,'defaultTextFontName', 'Garamond')
set(groot,'defaultAxesTickLabelInterpreter','latex'); 
set(groot,'defaultLegendInterpreter','latex');
cd(pwd); 
addpath(strcat(pwd,'/auxiliary_functions/'));
project_folder                          = fileparts(fileparts(pwd));
root_folder                             = fileparts(project_folder);

%%%%%%%%%%%%%%%
% Model setup
%%%%%%%%%%%%%%%

% Structural parameters
pars.r                                  = -1/12*log(0.70);                                                              % Monthly model with annual discount rate of 0.70 
pars.k                                  = 0.163;                                                                        % Upper bound on adjustment intensity
pars.M_c                                = 1;                                                                            % Cash stock average
pars.M_e                                = 0.97;                                                                         % Electronic money    
pars.C                                  = 0.060;                                                                        % Strength of complementarities           
pars.mean_reversion                     = 120;                                                                          % Shock persistence, expressed as the number of days to get back to within 10% of initial value
pars.theta                              = -log(1-0.90)/( pars.mean_reversion/30 );                                      % Shock persistence, expessed as OU parameter
pars.sigma                              = 0.060;                                                                        % Standard deviation of money shock    
pars.T                                  = 12*10;                                                                        % Date after which mean-reversion goes to zero
pars.tIRF                               = 12*2;                                                                         % IRF horizon

% Discretization parameters
Nu.NX                                   = 1e2+1;                                                                        % # of gridpoints for X
Nu.N                                    = 1e2;                                                                          % (1/2)*( # of gridpoints for M - 1 )
Nu.Nmin                                 = ceil(pars.theta/pars.sigma^2* ...
                                               ((pars.r+pars.k+pars.theta)/(pars.r+pars.k))^2* ...
                                              max((pars.M_c-pars.M_e)^2,(pars.M_e + pars.C - pars.M_c)^2)); 
Nu.N                                    = max(Nu.N,Nu.Nmin); 
Nu.Deltat                               = 1/(Nu.N*pars.theta);                                                          % Time period, expressed as fraction of a month
Nu.theta_target                         = pars.theta;
Nu.solutionmode                         = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% C > 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pars.C                                  = 0.060;
solution                                = solve(pars,Nu,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% C = 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pars.C                                  = 0;
solution_noC                            = solve(pars,Nu,1);

%%%%%%%%%%
% Plots
%%%%%%%%%%
fignum                                  = 1;
tvec                                    = (-solution.preT.Nu.Deltat):solution.preT.Nu.Deltat:solution.preT.pars.T;
tvecIRF                                 = (-solution.preT.Nu.Deltat):solution.preT.Nu.Deltat:solution.preT.pars.T;
color1                                  = [0, 0.4470, 0.7410];
color2                                  = [0.6350, 0.0780, 0.1840];
gray                                    = [0,0,0]+0.3;
gray2                                   = [0,0,0]+0.5;
gray3                                   = [0,0,0]+0.7;
color3                                  = [0.9290, 0.6940, 0.1250];
color4                                  = [0.4660, 0.6740, 0.1880];

% Stationary distribution - C = 0
ff1 = figure(fignum);
fignum = fignum+1;

    p2=area(solution_noC.preT.Nu.X,sum(solution_noC.preT.mu,1),'FaceColor',color1,'EdgeColor',color1,'FaceAlpha',0.5,'EdgeAlpha',1); 
    grid on;
    xlabel('$X_t$','interpreter','latex');
    %title('Stationary distribution of $X_t$','interpreter','latex','FontSize',12);
    ylim([0,0.12]);
    set(findall(gcf,'-property','FontSize'),'FontSize',20) 

    % Print   
    set(ff1,'Units','inches');
    screenposition                              = get(ff1,'Position');
    set(ff1,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff1.PaperPositionMode                       ='auto';
    exportgraphics(ff1,[root_folder '/output/figures/' 'Figure_2_2.pdf'],'BackgroundColor','none') 

% Phase diagram - C = 0 - Large shock
ff2 = figure(fignum);
fignum = fignum+1;

    X                                   = solution_noC.preT.Nu.X';
    Phi                                 = solution_noC.preT.Phi_0;
    Phiref                              = Phi(1);
    ylim1                               = 0.8;
    ylim2                               = 1.2;
    Deltat                              = solution_noC.preT.Nu.Deltat;
    k                                   = solution_noC.preT.pars.k;
    theta                               = solution_noC.preT.pars.theta;
    M_c                                 = solution_noC.preT.pars.M_c;
    NX_q                                = 1e1;
    NM_q                                = 1e1;
    Xgrid                               = repmat(linspace(0,1,NX_q),[NM_q 1]);
    Mgrid                               = repmat(linspace(ylim1,ylim2,NM_q)',[1 NX_q]);
    Phi                                 = solution_noC.preT.Phi_0;
    Phigrid                             = interp1(X,Phi,Xgrid(:));
    Phigrid                             = reshape(Phigrid,size(Xgrid));
    dXgrid                              = k*Deltat*( (1-Xgrid).*( Mgrid <= Phigrid ) - Xgrid.*( Mgrid > Phigrid ) );
    dMgrid                              = theta*Deltat*( M_c - Mgrid );
    index_above                         = Mgrid > Phigrid;
    index_below                         = Mgrid <= Phigrid;

    X                                   = solution_noC.preT.Nu.X';
    Phi                                 = solution_noC.preT.Phi_0;
    Phiref                              = Phi(1);
    tPlot                               = solution.IRF.tvec;
    ylim1                               = 0.8;
    ylim2                               = 1.2;
    M_0                                 = pars.M_c;
    M_t                                 = pars.M_c*(1-exp(-pars.theta*tPlot)) + exp(-pars.theta*tPlot)*(1-0.175)*pars.M_c;
    X_0                                 = 0;
    X_t                                 = exp( - (pars.k)*tPlot )*X_0 + (1-exp( - pars.k*tPlot ));
    indexswitch                         = find(M_t >= Phiref,1,'first');
    X_t(indexswitch:end)                = exp( - (pars.k)*(tPlot(indexswitch:end)-tPlot(indexswitch)) )*X_t(indexswitch);

    
    p2=plot(X,Phi,'-','color',color1,'linewidth',0.5); 
    hold on;
    plot(X,pars.M_c*ones(size(Phi)),'--','color',gray3,'linewidth',0.5);
    hold on;
    ylim([ylim1,ylim2]);
    xlabel('$X_t$','interpreter','latex');
    ylabel({'$M_t \quad$ ',''},'interpreter','latex','Rotation',0);
    th1 = text(1+0.04,pars.M_c,{'$\quad M^c$'},'interpreter','latex','HorizontalAlignment','center','color',gray3,...
                                'VerticalAlignment','middle','fontsize',12); 
    th2 = text(1+0.09,Phi(end),{'$\mathbf{\Phi}_t^{(0)}$'},'color',color1,'interpreter','latex','HorizontalAlignment','center',...
                                'VerticalAlignment','middle','fontsize',12); 
    set(th2,'Rotation',0)
    ax = gca;
    ax.YLim = [ylim1 ylim2];
    ax.YTick = [];
    ax.XLim = [0.00 1.00];
    rna = fill([X' flip(X)'],[(Phi)' (ax.YLim(2)*ones(size(X)))'],color3,'LineStyle','none');
    rna.FaceAlpha = 0.15;
    rna = fill([X' flip(X)'],[(Phi)' (ax.YLim(1)*ones(size(X)))'],color4,'LineStyle','none');
    rna.FaceAlpha = 0.15;
    text(0.5,Phiref-0.05,{'$dX_t > 0$'},'color',color4,'interpreter','latex',...
                    'HorizontalAlignment','center','Fontsize',10);  

    text(0.5,Phiref+0.20,{'$dX_t < 0$'},'color',color3,'interpreter','latex',...
                    'HorizontalAlignment','center','Fontsize',10); 

    plot(X_0,M_0,'o','MarkerSize',10,'MarkerFaceColor','white','MarkerEdgeColor',color1); hold on;
    plot(X_t(1),M_t(1),'o','MarkerSize',10,'MarkerFaceColor',color1,'MarkerEdgeColor',color1); hold on;
    set(findall(gcf,'-property','FontSize'),'FontSize',25)    
    
    % Print
    set(ff2,'Units','inches');
    screenposition                              = get(ff2,'Position');
    set(ff2,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff2.PaperPositionMode                       ='auto';
    set(ff2,'position',[0,0,8,7])
    X_t                                         = X_t(1:(indexswitch-1));
    M_t_temp                                    = M_t(1:(indexswitch-1));
    p1 = plot(X_t,M_t_temp,'-','color',color1,'LineWidth',2,'MarkerSize',0.01,'MarkerFaceColor',gray2,'MarkerEdgeColor',color1); hold on;
    X_t                                         = exp( - (pars.k)*tPlot )*X_0 + (1-exp( - pars.k*tPlot ));
    X_t(indexswitch:end)                        = exp( - (pars.k)*(tPlot(indexswitch:end)-tPlot(indexswitch)) )*X_t(indexswitch);
    p2 = plot(X_t,M_t,'-','color',color1,'LineWidth',2,'MarkerSize',0.01,'MarkerFaceColor',gray2,'MarkerEdgeColor',color1); hold on;
    line2arrow(p1);
    line2arrow(p2);
    exportgraphics(ff2,[root_folder '/output/figures/' 'Figure_2_1.pdf'],'BackgroundColor','none')

% Stationary distribution - C > 0
ff3 = figure(fignum);
fignum = fignum+1;

    p2=area(solution.preT.Nu.X,sum(solution.preT.mu,1),'FaceColor',color2,'EdgeColor',color2,'FaceAlpha',0.5,'EdgeAlpha',1); 
    grid on;
    xlabel('$X_t$','interpreter','latex');

    % Print  
    set(ff3,'Units','inches');
    screenposition                              = get(ff3,'Position');
    set(ff3,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff3.PaperPositionMode                       ='auto';
    set(findall(ff3,'-property','FontSize'),'FontSize',20) 
    exportgraphics(ff3,[root_folder '/output/figures/' 'Figure_2_4.pdf'],'BackgroundColor','none') 

% Phase diagram - C > 0 - large shock
ff4 = figure(fignum);
fignum = fignum+1;

    X = solution.preT.Nu.X';
    Phi = solution.preT.Phi_0;
    ylim1 = 0.8;
    ylim2 = 1.2;

    X                                   = solution.preT.Nu.X';
    Phi                                 = solution.preT.Phi_0;
    Phiref                              = Phi(1);
    tPlot                               = solution.IRF.tvec;
    ylim1                               = 0.8;
    ylim2                               = 1.2;
    M_0                                 = pars.M_c;
    M_t                                 = pars.M_c*(1-exp(-pars.theta*tPlot)) + exp(-pars.theta*tPlot)*(1-0.175)*pars.M_c;
    X_0                                 = 0;
    X_t                                 = exp( - (pars.k)*tPlot )*X_0 + (1-exp( - pars.k*tPlot ));
    PhiX_t                              = interp1(X',Phi,X_t);
    indexswitch                         = find(M_t >= PhiX_t,1,'first');
    if ~isnan(indexswitch)
        X_t(indexswitch:end)                = exp( - (pars.k)*(tPlot(indexswitch:end)-tPlot(indexswitch)) )*X_t(indexswitch);
    end

    p2=plot(X,Phi,'-','color',color2,'linewidth',0.5); 
    hold on;
    plot(X,pars.M_c*ones(size(Phi)),'--','color',gray3,'linewidth',0.5);
    hold on;
    ylim([ylim1,ylim2]);
    xlabel('$X_t$','interpreter','latex');
    ylabel({'$M_t \quad$',''},'interpreter','latex','Rotation',0);
    th1 = text(1+0.07,pars.M_c,{'$M^c$'},'interpreter','latex','HorizontalAlignment','center','color',gray3,...
                                'VerticalAlignment','middle','fontsize',12); 
    th2 = text(1+0.12,Phi(end),{'$\mathbf{\Phi}_t(X_t)$'},'color',color2,'interpreter','latex','HorizontalAlignment','center',...
                                'VerticalAlignment','middle','fontsize',12); 
    set(th2,'Rotation',0)

    ax = gca;
    ax.YLim = [ylim1 ylim2];
    ax.YTick = [];
    ax.XLim = [0.00 1.00];
    
    rna = fill([X' flip(X)'],[(Phi)' (ax.YLim(2)*ones(size(X)))'],color3,'LineStyle','none');
    rna.FaceAlpha = 0.15;

    rna = fill([X' flip(X)'],[(Phi)' (ax.YLim(1)*ones(size(X)))'],color4,'LineStyle','none');
    rna.FaceAlpha = 0.15;

    text(0.5,Phiref-0.05,{'$dX_t > 0$'},'color',color4,'interpreter','latex',...
                    'HorizontalAlignment','center','Fontsize',10);  

    text(0.5,Phiref+0.20,{'$dX_t < 0$'},'color',color3,'interpreter','latex',...
                    'HorizontalAlignment','center','Fontsize',10);  

    plot(X_0,M_0,'o','MarkerSize',10,'MarkerFaceColor','white','MarkerEdgeColor',color2); hold on;
    %text(X_t(1)-0.12,M_t(1)-0.01,'$M_0$','color','white','interpreter','latex','fontsize',14); hold on;
    plot(X_t(1),M_t(1),'o','MarkerSize',10,'MarkerFaceColor',color2,'MarkerEdgeColor',color2); hold on;
    %text(X_t(1)-0.12,M_t(1)-0.01,'$M_0$','color',color2,'interpreter','latex','fontsize',14); hold on;
    set(findall(ff4,'-property','FontSize'),'FontSize',25)    
    set(ff4,'Units','inches');
    screenposition                              = get(ff4,'Position');
    set(ff4,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff4.PaperPositionMode                       = 'auto';
    set(ff4,'position',[0,0,8,7])
    p4=plot(X_t,M_t,'->','color',color2,'LineWidth',1,'MarkerSize',0.01,'MarkerFaceColor',gray2,'MarkerEdgeColor',color2); hold on;
    line2arrow(p4);
    exportgraphics(ff4,[root_folder '/output/figures/' 'Figure_2_3.pdf'],'BackgroundColor','none');  

